##Authors and contact: Naielly Lopes Marques (naielly.lopes@iag.puc-rio.br), Carlos de Lamare Bastian-Pinto (carbastian@gmail.com), and Luiz Eduardo Teixeira Brand�o (brandao@iag.puc-rio.br). Institution: Pontifical Catholic University of Rio de Janeiro (PUC-RJ). Department: IAG Business School
##Purpose: This R code helps researchers and practitioners model project cash flows for real option applications considering the correct volatility estimation (Brand�o et al., 2012), dividend yield modeling and lattice building
##Link to published paper: A Tutorial for Modeling Real Options Lattices from Project Cash Flows
##Link to code: https://doi.org/10.5281/zenodo.3885925
##Last update: June 9th, 2020

##Package needed for this code
if (!require(DescTools)) { 
  install.packages("DescTools") #We will use the NPV function of this package to calculate the Net Present Values
}
if (!require(fOptions)) { 
  install.packages("fOptions") #We will use the BinomialTreePlot function of this package to plot the lattices
}

##Modelling price
#Parameters - Here, you can change the input values to suit your project
n <- 10 #Depreciation duration time
a <- 0.03 #Price growth rate 
vol <- 0.15 #Price volatility 
P0 <- 100 #Price at t = 0
nt <- 10000 #Number of simulations
#Calculations derived from the input values
i <- n #Number of time intervals
dt <- n/i #Time interval
t <- seq(from=dt,to=n,by=dt)
l <- length(t)
at <- a*dt
volt <- vol*sqrt(dt)
as <- at-(volt^2)/2
#Simulating price
#Creating a function to simulate a GBM process
GBM_function <- function(start, nsim, n, drift, volatility, growth) {
  x = matrix(NA, nrow=nsim, ncol=(n+1))
  x[,1] <- start
  for (j in 1:nsim) {
    x[j,2] <- start*exp(drift+volatility*rnorm(1))
  }
  for(j in 1:nsim){
    for (i in 3:(n+1)) {
      x[j,i] <- x[j,(i-1)]*exp(growth) 
    }
  }
  return(x)
}
X <- GBM_function(P0,nt,l,as,volt,at)
X #This will generate matrix X of dimension nt x (l+1), where the first column represents the prices at t = 0 (P0 = 100) and the other columns the prices simulated at each time point until the tenth year of the project

##Calculating Cash Flow
#Parameters - Here, you can change the input values to suit your project
r <- 0.06 #Risk free rate
k <- 0.12 #Discount rate
g <- 0.03 #Perpetuity growth rate
prod <- 10000 #Production
VC <- 0.55 #Variable costs
FC <- 300000 #Fixed costs
I <- 1500000 #Investment
EI <- 50000 #Extra investments
IT <- 0.34 #Income tax
#Calculations derived from input values
rt <- r*dt
kt <- k*dt
gt <- g*dt
R <- prod*X[,-1] #Revenue
FC <- matrix(rep(FC),nrow=nt,ncol=l) #Fixed costs matrix
RO <- R-R*VC-FC #Operating revenue
Imatrix <- matrix(rep(I),nrow=nt,ncol=l) #Investment matrix
EImatrix <- matrix(rep(EI),nrow=nt,ncol=l) #Extra investments matrix
Dep0 <- Imatrix[,1]/l #Depreciation at t = 1
Dep <- matrix(rep(0),nrow=nt,ncol=(l-1)) #Depreciation
for (i in 1:(n-1)) {
  Dep[,i] <- (Imatrix[,i]+EImatrix[,i]*i)/l
}
Dep <- cbind(Dep0,Dep)
colnames(Dep) <- NULL
EBIT <- RO-Dep #Earnings Before Interest and Taxes
#Cash Flow Matrix
FCF <- EBIT-IT*EBIT-EImatrix+Dep #Free Cash Flow
Perp <- FCF[,l]/(kt-gt)*(1+gt)*(1+kt) #Perpetuity
FCF <- cbind(FCF,Perp)
colnames(FCF) <- NULL #This will generate a cash flow matrix (FCF) of dimensions nt x l. In addition, since we consider that this project has a continuation value (perpetuity), we have included a column in this matrix that represents the perpetual cash flows of that project

##Calculating PV ex ante and ex post
PV <- matrix(rep(0),nrow=nt,ncol=1) #Present Value
for (i in 1:nt) {
  PV[i,] <- NPV(kt,FCF[i,],seq(along=FCF[i,]))
}
NPV <- PV-Imatrix[,1] #Net Present Value
PVa <- matrix(rep(0),nrow=nt,ncol=l) #Ex ante Present Value
for(j in 1:nt){
  for (i in 1:l){
    PVa[j,i] <- NPV(kt,FCF[j,((i+1):(l+1))],seq(along=FCF[j,((i+1):(l+1))]))+FCF[j,i]
    #This will generate a matrix of Ex ante Present Values
  }
}
PVp <- matrix(rep(0),nrow=nt,ncol=l) #Ex post Present Value
for (i in 1:l) {
  for(j in 1:nt){
    PVp[j,i] <- NPV(kt,FCF[j,((i+1):(l+1))],seq(along=FCF[j,((i+1):(l+1))]))
    #This will generate a matrix of Ex post Present Values
  }
}

##Calculating the project volatility
PVd <- mean(PV)
PVd #Here, we find that the average project value is: V0 = $1,661,448, yielding a Net Present Value (NPV) of $161,448
#This value can vary slightly as it is the result of a Monte Carlo Simulation
Ret <- (PVa[,1]/PVd) #Return
lRet <- log(Ret)
sig <- sd(lRet,na.rm = TRUE)
sig #We find that the project volatility is 33%
#This value can vary slightly as it is the result of a Monte Carlo Simulation

##Ex ante and ex post Lattices
u <- exp(sig) #Upside multiplying factor
d <- 1/u #Downside multiplying factor
p <- (1+rt-d)/(u-d) #Probability
divr <- FCF[,-(l+1)]/PVa #Dividend rate
div <- cbind(1,(1-divr)) #Dividends
Latta <- matrix(NaN,(l+1),(l+1)) #Ex ante Lattice, see Apendix I of the published paper
Lattp <- matrix(NaN,(l+1),(l+1)) #Ex post Lattice, see Apendix I of the published paper
Lattp[1,1] <- PVd
Latta[1,2] <- PVd*u
for (a in 2:(l+1)) {
  for (b in 3:(l+1)) {
    Lattp[1,a] <- Latta[1,a]*div[1,a]
    Latta[1,b] <- Lattp[1,(b-1)]*u
  }
}
for (b in 2:(l+1)) {
  for (a in 2:(l+1)) {
    Latta[b,a] <- Lattp[(b-1),(a-1)]*d
    Lattp[b,a] <- Latta[b,a]*div[1,a]
  }
}
#Saving ex ante and ex post lattices in a pdf file 
pdf(file = "V Lattice with ex ante values.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
BinomialTreePlot(Latta, dy = 1, cex = 0.8, axes = FALSE, ylim = c(-15, 15),
                 xlab = " ", ylab = " ") #We use this function of fOptions package to plot the ex ante lattice
title(main = "V Lattice with ex ante values")
dev.off()
pdf(file = "V Lattice with ex post values.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
BinomialTreePlot(Lattp, dy = 1, cex = 0.8, axes = FALSE, ylim = c(-15, 15),
                 xlab = " ", ylab = " ") #We use this function of fOptions package to plot the ex post lattice
title(main = "V Lattice with ex post values")
dev.off()

##Incorporating abandonment and expansion options
#To find the residual value of the project in case of abandonment, we calculate the depreciated asset value in each year until year n by discounting the depreciation amount from the investments and we multiply these values by the abandonment factor
Dep0 <- c(I/l) #Depreciation at t = 1
Dep <- matrix(rep(0,(l-1)),nrow=1) #Depreciation
for (i in 1:(l-1)) {
  Dep[,i] <- (I+EI*i)/l
}
Dep <- c(Dep0,Dep)
EIvector <- (rep(EI,l))
Depasset <- (rep(NaN,(l-1))) #Depreciated asset
for (i in 2:n) {
  Depasset[(i-1):(l-1)] <- (I+sum(EIvector[-c(i:n)])-sum(Dep[-c(i:n)]))
}
Depasset <- cbind(0, t(as.matrix(Depasset, ncol = (l-1), nrow = 1)), (I + sum(EIvector[1:n]) - sum(Dep[1:n])))
#Parameters - Here, you can change the input values to suit your project
abandf <- 0.8 #Abandonment factor 
expf <- 1.8 #Expansion factor 
expc <- 1200000 #Expansion cost
#Calculation derived from input values
Residual <- Depasset*(abandf) #Residual value

##Ex ante and ex post Lattices with options
#Here, the code evaluates backwards the maximum value between maintaining, abandoning and expanding the project each year until year n
Lattpo <- matrix(NaN,(l+1),(l+1)) #Ex post Lattice with Options, see Apendix I of the published paper
for (i in 1:(l+1)) {
  Lattpo[i,(l+1)] <- max(Lattp[i,(l+1)],Lattp[i,(l+1)]*expf-expc,Residual[,(l+1)])
}
Lattao <- matrix(NaN,(l+1),(l+1)) #Ex ante Lattice with Options, see Apendix I of the published paper
for (i in 1:(l+1)) {
  Lattao[i,(l+1)] <- Lattpo[i,(l+1)]+(Latta[i,(l+1)]-Lattp[i,(l+1)])
}
for (j in l:1) {
  for (i in 1:l) {
    Lattpo[i,j] <- max(Lattp[i,j]*expf-expc,Residual[,j],(Lattao[i,(j+1)]*p+Lattao[(i+1),(j+1)]*(1-p))/(1+rt))
    Lattao[i,j] <- Lattpo[i,j]+(Latta[i,j]-Lattp[i,j])
  }
}
Lattpo[1,1] #At the starting step of this lattice (Lattpo), we can verify that the project value considering the expansion and abandonment options is: $2,109,671
#This value can vary slightly as it is the result of a Monte Carlo Simulation
#Saving ex ante and ex post lattices with exercise of real options in a pdf file 
pdf(file = "V Lattice with ex ante values and exercise of real options.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
BinomialTreePlot(Lattao, dy = 1, cex = 0.8, axes = FALSE, ylim = c(-15, 15),
                 xlab = " ", ylab = " ") #We use this function of fOptions package to plot the ex ante lattice with options
title(main = "V Lattice with ex ante values and exercise of real options")
dev.off()

pdf(file = "V Lattice with ex post values and exercise of real options.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
BinomialTreePlot(Lattpo, dy = 1, cex = 0.8, axes = FALSE, ylim = c(-15, 15),
                 xlab = " ", ylab = " ") #We use this function of fOptions package to plot the ex post lattice with options
title(main = "V Lattice with ex post values and exercise of real options")
dev.off()

##Graphs of ex post lattices with and without options
#First, we create a function to mirror any lattice so that we can plot it
latt_inverter <- function(vec){
  vec2 <- vec[!is.na(vec)]
  vec3 <- c(rep(NaN,length(vec)-length(vec2)),vec2)
  vec3
} 
Lattpb <- apply(Lattp,2,latt_inverter) #Mirrored Ex post Lattice
x <- c(0:l) #x-axis corresponds to the project's duration years
pdf(file = "Ex post lattices with and without abandonment and expansion options.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
par(mar=c(5.5,5.5,5.5,5.5)) #You can choose the plot margins
for (b in 1:(l+1)) {
  y <- (Lattp[b,]/1000) #y-axis corresponds to the ex post lattice without options
  plot(x,y,type="l",ylim=c(min(Lattp/1000,na.rm=TRUE),max(Lattpo/1000,na.rm=TRUE)),las=1,col="blue",ann = FALSE)
  par(new=T)
}
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpb[b,]/1000) #y-axis corresponds to the mirrored ex post lattice without options
  plot(x,y,type="l",ylim=c(min(Lattp/1000,na.rm=TRUE),max(Lattpo/1000,na.rm=TRUE)),las=1,col="blue",ann = FALSE)
}
par(new=T)
value <- round(y[1],digits=1)
text(0,rep(value-1000,1),(value),col=4)
Lattpob <- apply(Lattpo,2,latt_inverter) #Mirrored Ex post Lattice with Options
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpo[b,]/1000) #y-axis corresponds to the ex post lattice with options
  plot(x,y,type="l",ylim=c(min(Lattp/1000,na.rm=TRUE),max(Lattpo/1000,na.rm=TRUE)),las=1,col="red",ann = FALSE)
}
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpob[b,]/1000) #y-axis corresponds to the mirrored ex post lattice with options
  plot(x,y,type="l",ylim=c(min(Lattp/1000,na.rm=TRUE),max(Lattpo/1000,na.rm=TRUE)),las=1,col="red",ann = FALSE)
}
par(new=T)
value <- round(y[1],digits=1)
text(0,rep(value+1000,1),(value),col=2)
mtext(side = 1,text = "Year",line = 3) #x-axis label
mtext(side = 2, text = "Project Value ($)",line = 4) #y-axis label
dev.off()

##Plotting ex post lattice with and without options in log scale
pdf(file = "Ex post lattices with and without abandonment and expansion options in log scale.pdf", #You can name the file and change the directory you want to save it
    width = 13.00, #You can choose the width of the plot in inches
    height = 10.00) #and the height of the plot in inches
par(mar=c(5.5,5.5,5.5,5.5)) #You can choose the plot margins
marks <- c(10,100,1000,10000,100000)
for (b in 1:(l+1)) {
  y <- (Lattp[b,]/1000) #y-axis corresponds to the ex post lattice without options
  plot(x,y,type="l",log="y",ylim=c(10,100000),yaxt="n",col="blue",ann=FALSE)
  axis(2,at=marks,labels=format(marks,scientific=FALSE),las=2)
  par(new=T)
}
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpb[b,]/1000) #y-axis corresponds to the mirrored ex post lattice without options
  plot(x,y,type="l",log="y",ylim=c(10,100000),yaxt="n",col="blue",ann = FALSE)
  axis(2,at=marks,labels=format(marks,scientific=FALSE),las=2)
}
par(new=T)
value <- round(y[1],digits=1)
text(0,rep(value-500,1),(value),col=4)
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpo[b,]/1000) #y-axis corresponds to the ex post lattice with options
  plot(x,y,log="y",type="l",ylim=c(10,100000),yaxt="n",col="red",ann = FALSE)
  axis(2,at=marks,labels=format(marks,scientific=FALSE),las=2)
}
for (b in 1:(l+1)) {
  par(new=T)
  y <- (Lattpob[b,]/1000) #y-axis corresponds to the mirrored ex post lattice with options
  plot(x,y,log="y",type="l",ylim=c(10,100000),yaxt="n",col="red",ann = FALSE)
  axis(2,at=marks,labels=format(marks,scientific=FALSE),las=2)
}
par(new=T)
value <- round(y[1],digits=1)
text(0,rep(value+500,1),(value),col=2)
mtext(side = 1,text = "Year",line = 3) #x-axis label
mtext(side = 2, text = "Project Value ($)",line = 4) #y-axis label
dev.off()
